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Equilibrium Surface Current and Role of U(l) Symmetry: 
sum rule and snrface pertnrbations 
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We discuss effects of surface perturbations on equilibrium surface currents which contribute to 
orbital magnetization and orbital angular momentum in systems without time reversal symmetry. 
We show that, in a U(l) particle number conserving system, disorder and other perturbations at a 
surface do not affect the equilibrium surface current and corresponding orbital magnetization due to 
a sum rule which is analogous to Luttinger’s theorem. On the other hand, for a superfluid, the sum 
rule is no longer applicable and hence the surface mass current and corresponding orbital angular 
momentum can depend on details of a surface. 

PACS numbers: Valid PACS appear here 


I. INTRODUCTION 

The orbital magnetization (OM) and orbital angular 
momentum (0AM) are one of the most fundamental 
physical quantities in condensed matter physics, which 
arise in systems without time reversal symmetry due to 
equilibrium circulating currents. In a continuum (non¬ 
lattice) system which does not break translational sym¬ 
metry in bulk regions, an equilibrium circulating current 
flows only near its boundary. While in a lattice sys¬ 
tem, a circulating current is also possible around each 
atom in addition to a surface current. In spite of their 
physical importance, however, OM/OAM and associated 
equilibrium surface currents are not well understood. In¬ 
deed, the OM arising from surface currents are usually 
not taken into account but only locally circulating cur¬ 
rents around atoms are included, when one calculates 
total magnetization of a system in the density functional 
theoryi. This would be partly because it is difficult to ap¬ 
propriately treat surface currents theoretically and also 
effects of surface currents on the total magnetization are 
expected to be small in conventional ferromagnets such 
as Fe and Ni. 

However, contributions of surface currents to 
OM/OAM would become important in several interesting 
systems such as topological insulators/superconductors 
without time reversal symmetry and fractional quantum 
Hall systems where there are chiral edge modes^Ti^. In 
the last decade, many theoretical studies for OM have 
been proposed for band insulators and metals including 
topological systems, and they give a beautiful formula 
for calculating OMi^“— . This formula involves only the 
Bloch wavefunctions which are bulk properties of a sys¬ 
tem, and hence, the OM and associated surface currents 
can be regarded as bulk quantities. These results are 
in good agreement with our general expectations that 
magnetism is a bulk property which is independent of 
surface details in real materials. Conceptually, they 
can be thought as a nice realization of the bulk-surface 
correspondence in a general sense that bulk properties 
are determined by surface physics and vice versa. At 
the same time, however, one may naively expect that 


surface currents could be affected by perturbations near 
surfaces which should exist in real materials, such as 
surface disorder, deformation of Wannier functions, local 
inversion symmetry breaking, weak screening of the 
Coulomb interaction, and so on. 

A partial answer to this fundamental question can be 
obtained by following the derivations of the formula for 
OM. The formula has been derived in three different 
ways; (i) direct calculations of circulating currents for 
trivial band insulators in the presence of boundaries^i^, 

(ii) semi-classical wavepacket approximations — , and 

(iii) taking derivative of free energy with respect to mag¬ 
netic fields under the periodic boundary condition 
Following the derivation (i), one could find that the sur¬ 
face currents are independent of details of the boundaries 
when the ground state is a simple band insulator. Based 
on this, it was argued that effects of surface perturba¬ 
tions on OM in such simple insulators are irrelevant^. 
Although the discussion presented in the derivation (i) 
cannot be applied to other systems, the same formula 
was obtained by the derivations (ii) and (iii). The sur¬ 
face currents are found to be independent of gradient of 
surface potentials within the semi-classical wavepacket 
approximations in the derivation (ii). In this approxi¬ 
mation, however, the length scale of surface potentials 
should be much longer than the wavepacket size. In the 
derivation (iii), although the OM is given as a bulk quan¬ 
tity by its definition for systems with periodic boundary 
conditions, connections to a surface current which exists 
in a realistic finite size system are unclear— There¬ 
fore, in spite of the surprisingly convincing agreements 
between the three derivations, it is still not clear why the 
OM and corresponding surface currents are given as bulk 
quantities which are independent of surface details. 

In contrast to non-superconducting systems discussed 
above, surface perturbations do become relevant for 
0AM in chiral superfluids which break time reversal sym¬ 
metry. The 0AM in chiral superfluids has been a long¬ 
standing issue and attracting much interest since the dis¬ 
covery of ^He A-phase— . Effects of surface roughness 
have been studied in the context of ^He A-phase and 
Sr 2 Ru 04 and it was theoretically argued that the 


surface mass currents and resulting 0AM in chiral su¬ 
perfluids depend on surface roughness in weak coupling 
region o^^i^^'^^ . In case of domain walls between chiral 
BCS superfluids with opposite chiralities, the boundary 
currents strongly depend on difference in U(l) phases 
of order parameter o^^'^^ . Besides, when the surface is 
sharp, net surface currents vanish for higher order pair¬ 
ing states such as d+id,f + if-wave BCS superfluids 
due to hidden depairing effects which exist even for clean 
surfaces^i^. These results suggest that, in contrast to 
U(l) symmetric systems, surface currents and 0AM in 
superfluids with broken time reversal symmetry depend 
on surface details and are not bulk quantities. However, 
physical origins of the reduction of the surface currents 
and difference between the superfluids and U(l) sym¬ 
metric systems have not been well understood. Besides, 
while surface perturbations are known to be relevant in 
the weak coupling BCS states where Cooper pairs are ex¬ 
tended in space, they have not been discussed so far for 
the strong coupling BEC states where Cooper pairs are 
strongly bounded^. For such tightly bounded pairing 
states, one may naively expect that the surface currents 
are robust against surface perturbations. 

In this paper, in order to clarify the different behav¬ 
iors of the surface currents and corresponding OM/OAM 
in systems with or without U(l) symmetry, we discuss 
effects of surface perturbations from a general point of 
view. Based on a sum rule argument which is analogous 
to Luttinger’s theorem— and numerical calculations, 
we show that surface currents are robust against sur¬ 
face perturbations in general U(l) symmetric systems. 
On the other hand, for superfluids without time reversal 
symmetry, the sum rule is no longer applicable and sur¬ 
face currents can depend on surface perturbations such 
as surface roughness. Especially, it is shown that surface 
currents are suppressed in chiral superfluids on a lattice 
even for strong coupling BEC states. 

This paper organizes as follows. In Sec. m we discuss 
surface currents in U(l) symmetric systems. We show a 
sum rule for the surface current in general systems and 
confirm it by numerical calculations for a concrete model. 
The OM formula is revisited based on the sum rule argu¬ 
ment. Then, surface currents in U(I) broken systems are 
examined in Sec lIIIl Finally, we summarize this paper in 
Sec. IS 


II. SYSTEM WITH U(l) SYMMETRY 
A. sum rule argument 

In this section, we discuss surface currents in U(I) sym¬ 
metric systems. For simplicity, we consider 2-dimensional 
lattice models whose size is x Ny with the open bound¬ 
ary condition for a;-direction and the periodic boundary 
condition for j/-direction as shown in Fig. [T] Our argu¬ 
ment holds also for continuum models as discussed in 
Appendix A. We assume that there is no time reversal 
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FIG. 1: Schematic picture of the two-dimensional cylinder 
system. 


symmetry and there can exist equilibrium surface cur¬ 
rents for the models considered. Our Hamiltonian reads 
in general 

H = ^ ^ (fcy (fcy ) -f Hint T Hsurt (1) 

where Cxiiky) is an annihilation operator of electrons at 
position X, with wavenumber along y-direction ky and 
other quantum number I such as orbitals and spins. 
The matrix K describes one-particle Hamiltonian. We 
assume that the interaction H^t is simply short-range 
density-density interactions or on-site interactions so that 
it commutes with the total particle density operator, 
[ui^Hint] = 0. Hsurf is the surface perturbation term 
which is finite only near the surface and zero otherwise, 
and is also assumed to satisfy [ui, Hsiuf] = 0. Therefore, 
the U(l) current operator is simply given by the kinetic 
term only through the continuity equation in the present 
study. 

In this geometry, the equilibrium surface current at the 
left (right) surface is given by 

ly^ ^ ~ '^V,xl,x'l'(ky){(^^i{ky)Cx'l'{ky))^ 

y xx'eSLiSR) 

( 2 ) 

in unit of the electron charge e. is a region only 

near the left (right) surface whose width is much smaller 
than the system width Nx- Since effects of a surface 
should not propagate deep into bulk regions, we can de¬ 
fine such regions in general systems^. Vy is a velocity 
matrix given by Vy(ky) = dK{ky)/dky. In the absence of 
Hsurf it is obvious that Iy°^ = ly 1:^ = 0 under a nat¬ 
ural assumption that the Hamiltonian m has inversion 
symmetry in the x-direction, x •<->■ —x. Since current den¬ 
sity which contributes to the surface current is localized 
only around the surface and it vanishes in bulk regions 
out of 5 'l the total surface current is written as 

ly = Vy;xl,x'l' iky){(^^i{ky)Cx'l' {ky)) 

^ xx'eall sites 

= iuj)e^^°^]. (3) 
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Here, G{k) = G{ky,iuj) = —{{c{ky)c^ {ky))){iijj) is the 
matrix Matsubara Green’s function, and trace describes 
summation over all the indices, x, I and w. We note 
that, although off diagonal elements become finite if the 
system breaks some symmetries together with time re¬ 
versal symmetry such as spin rotation symmetry in a 
ferromagnetic state, our discussion holds in a parallel 
way by introducing order parameters into G and mod¬ 
ifying the definition of Vy appropriately. For example, if 
the translational symmetry along the y-direction is pre¬ 
served, only possible off diagonal elements in the Green’s 
function are of the form {cx'i'{ky)cl.i{ky)) which have 
already been included in the above expression. If the 
translational symmetry along the y-direction is broken, 
although off diagonal elements with different momenta 
ky and ky + Qy with an ordering vector Qy become 
finite, Eq. m is still valid after modifying Vy{ky) to 
Vy{ky) = disLg{vy{ky),Vy{ky + Qy)) whcrB ky bolongs to 
the reduced Brillouin zone. 

Now we turn on i/gurf which is finite only near the left 
surface and zero in all the other regions including the 
right surface. Our fundamental assumption is that is 
unchanged by the left surface perturbation Ffgurf as long 
as the width of the cylinder is large enough, which means 
that if changes, it is totally attributed to change in 
ly . We also assume that i?surf is translationally symmet¬ 
ric along the surface. In case of surface disorder, ifgurf is 
still translationally symmetric on average, for which the 
disorder-averaged Green’s function is diagonal with re¬ 
spect to ky. Then, all the effects of such perturbations in 
one-particle quantities can be incorporated into the (av¬ 
eraged) selfenergy matrix E(fc). Effects of iJint are also 
described by the selfenergy, and if the system is in bro¬ 
ken symmetry states, we include the corresponding order 
parameter matrix in E(/c). The Green’s function G{k) is 
now written as a matrix inverse of [iuj — K{ky) — E(fc)]. 
Therefore, 
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( 4 ) 


It is clear that the first term vanishes. (If the transla¬ 
tional symmetry is broken, ky = ±7r should be replaced 
by appropriate boundary fc-vectors of the reduced Bril¬ 
louin zone.) Second term also vanishes because of the 
Luttinger-Ward identit y"^^ ^ 
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= 0 . 


( 5 ) 


Although the Luttinger-Ward identity may be violated in 
Mott insulators ^^4^ , it holds in general gapless/gapped 
Fermi liquids where there is no zero (pole) in the Green’s 
function (selfenergy), and also in weakly disordered sys¬ 
tems^. As long as the Luttinger-Ward identity is sat¬ 
isfied, -h /« = 0 even in the presence of the 

perturbations around the left surface, which means that 
ly is unchanged although the current density at each site 
could be affected by Hsmf- 

There is a nice analogy between the conservation of 
surface currents in the presence of surface perturbations 
and the well known conservation law, the Luttinger’s the- 
orerr^i"— . The Luttinger’s theorem claims that, though 
shape of a Fermi surface can be changed by interactions, 
its total volume is unchanged. Here, although a real 
space profile of U(l) surface current density ji can be 
modified by the surface perturbations, its sum around the 
surface I = X^iGsurf unchanged. In this view point, 
the conservation of the U(l) surface current in the pres¬ 
ence of surface perturbations is understood as a result 
of non-trivial cancellations of changes in ji at each site. 
This sum rule argues that surface currents and corre¬ 
sponding OM are bulk quantities which are independent 
of surface details as suggested in the modern theories on 
OMi^— , and this supports existence of the bulk-surface 
correspondence for these quantities. We note that, the 
sum rule is useful not only for developments in under¬ 
standing fundamental aspects of surface currents, but 
also for practical calculations of them. Although realistic 
surface potentials would be complicated in general, one 
can use a particular surface potential which is suitable 
for computing them, such as a hard wall potential and a 
sufficiently smooth potential. The sum rule guarantees 
that the calculation results of the surface current and OM 
are independent of the potential used. Indeed, the known 
formula for OM can be reproduced by calculating surface 
current contributions and also contributions from locally 
circulating currents under a sufficiently smooth surface 
potential, as will be discussed in Sec. Ill PI 

Up to now, we have not considered external or spon¬ 
taneously generated magnetic fields. In the presence 
of an applied magnetic field parallel to the z-direction, 
the hopping integrals acquire phase factors correspond¬ 
ing to the flux configuration by the Peierls substitution, 
and the Brillouin zone is reduced to a magnetic Bril¬ 
louin zone. Even in this case, by appropriately modi¬ 
fying the definition of v{ky), we can still derive a sum 
rule for a surface current in the same way. The surface 
current under a magnetic field for a time reversal sym¬ 
metric system is related to the Landau diamagnetism. It 
has been known that “bulk approaches” and “surface ap¬ 
proaches” are equivalent for the Landau diamagnetism; it 
can be calculated either by derivative of free energy with 
respect to the magnetic field under periodic boundary 
conditional^— or by computing a surface current in the 
presence of boundaries and these two approaches 

consistently give the same results. Experimentally, the 
skipping orbits which would be responsible for the sur- 
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face current have been observed in surface impedance 
measurements for many metals Our sum rule ar¬ 

gument would give a new understanding on the known 
equivalence between the two theoretical approaches. 

Finally, it is noted that, the sum rule can hold since 
the U(l) charges of the electrons are the well defined 
unique value e in the systems with U(l) symmetry. If 
U(l) symmetry is absent, however, particle sectors with 
charge +e and hole sectors with —e are mixed and the 
surface currents will not generally be conserved, as will 
be discussed in Sec. IIIII 


B. Bloch-Bohm’s theorem 

In this section, we discuss relations of our sum rule to 
the Bloch-Bohm’s theorem. The theorem states that net 
current should vanish in the ground states^^— . Here, 
we reexamine this theorem and point out that (i) it does 
not hold when surface currents are concerned and (ii) it 
needs some modifications when spontaneous symmetry 
breaking is involved even for currents running in bulk 
regions. 

We consider a general Hamiltonian as in Eq. o with¬ 
out symmetry breaking fields in a finite size cylinder LxL 
where open (periodic) boundary condition is imposed for 
x(y)-direction, and denote the ground state wavefunction 
as |0l). Since the system size is finite, there is no spon¬ 
taneous symmetry breaking and the ground state wave- 
function preserves the symmetries of the Hamiltonian in¬ 
cluding the time reversal symmetry if it is contained in 
H. When the Hamiltonian has time reversal symmetry, 
it is trivial that the total current vanishes because it 
is odd under the time reversal symmetry. In the follow¬ 
ing, we mainly discuss H with time reversal symmetry 
and H without time reversal symmetry will be touched 
on briefly. For time reversal symmetric H, we introduce 
external fields Ai/ex = A^ Aujiichcji' which break time 
reversal symmetry by assuming that the state |0l) of H 
has corresponding instability, where A is a dimensionless 
constant which should be taken as A —t- 0 in the end. 
Although we focus on U(l) symmetric external fields in 
this section, effects of U(I) symmetry breaking fields will 
be discussed in Sec. uni 

Following Bohm— , we consider a variational state 
= U0\Ol,\) where |0l,a) is the ground state of 
H\ = H 4-Ai/ex- The unitary operator is defined as Ug = 
exp[f0 and 9 = ^-KnjL, (n = 0, ±1, • • •) is 

required so that Ug is a well-defined operator. This re¬ 
quirement can easily be understood from the first quan¬ 
tization form of \Jg\ the twist operator for the wavefunc¬ 
tion 4'(ri, • • • , Tat) = '^{ri -\- L, - ■ ■ + L) is given 

by Ug = eyi\)[i9Y^-yj]. In order for .vm) = 

(C/6(4')(ri, • • • , r^r) to satisfy the periodic boundary con¬ 
dition, 0 = 27rn/T is needed. If 9 were not an integer 
multiple of the wavefunction (t/gik) is no longer 

an element of the domain of the Hamiltonian. This is 
essentially comes from the well-known fact that position 


operators cannot be well-defined on a torus, although it 
is sometimes missed even by experts^. Similarly in the 
second quantization formalism, c'j = UgCjUg = 
obeys the periodic boundary condition when 9 = 27rn/L 
is satisfied. Although this fundamental point has been 
missing in most of the previous studies concerning Bloch- 
Bohm’s theorem, this is important especially for dis¬ 
cussing possible net surface currents. 

We then evaluate energy difference between 10^, a) and 
|0l,a), 

5El,\ = {9l,\\H\\9l^x) — {0L,A|^fA|0L,A) 

+ \Y^{cos[9{yi - yg)] - l){0L,\\4Aa,jvCjv\QL,\) 

-f sin[6»(y, - yj)]{^L,\\c\iKajvCjv\{)L,\) 

+ A^isin[6»(yj - yj)\{GL,\\4iAiS'^3i'\^L,\)■ 

( 6 ) 

It is well known that the Lieb-Schultz-Mattis twist oper¬ 
ator Ug IS not helpful for higher dimensions other than 
one-dimension, in order to construct a variational state 
with low energy excitations. Indeed, since cos[0(?/i — 
vi)] — I ~ 0{\/L^) in the present model with a fi¬ 
nite hopping range It which is much shorter than L, 
the cos[9{yi — ?/j)]-terms are obviously 0(1) for large L 
limit. In general d-dimensional system of an isotropic 
size L for all directions, these terms are There¬ 

fore, the sin[d(?/i — ?/j)]-terms could become dominant 
for large L, if at least one of them is of order L or 
larger, namely ^zsin[0(2/i — yj)]{clKcj) ^ 0{L) or 
J2'isin[9{yi-yj)]{clAcj) 0{L). By expanding sin[-•• ], 
we see that one of the sin[- • • ]-terms is simply the total 
current running in the whole system, 0(Oi,.a| jyi\^L,\) 

in the leading order of 0 <C 1, where jyi is current den¬ 
sity. The original Bloch-Bohm’s theorem states that, 
since we can choose either 9 = 27r/L or 9 = —27r/L, 
in order for SEl,x to be non-negative, (Ol.a| 
must be zero in the thermodynamic limit L —>■ oo. How¬ 
ever, the Bloch-Bohm’s argument is based on an implicit 
assumption that the expectation value of the total cur¬ 
rent is 0{L‘^), and also the other sin-term 

is negligible, ^i{yi — yj)(c|Acj) ^ o(L). At the same 
tiem, we also should pay attention to the cos[- • • ]-terms 
which are 0(1). In the case of surface currents with¬ 
out any current density in the bulk, contributions only 
come from the surface regions Sl,r ~ 0(1), and there¬ 
fore 9{^-jyi) ^ (1/T) X L = 0(1) and it becomes the 
same order as the cos-terms. In d-dimensional systems, 
9{J2ijyi) ~ (1/T) X = 0{L‘^~‘^) and it competes 
with the cos-terms as well, for which we cannot immedi¬ 
ately conclude that finite (Ol,a| Jy*|Oi.A) 
contradicts with 5El^\. Hence, the Bloch-Bohm’s ar¬ 
gument is not applicable when surface currents are con¬ 
cerned, although it might give an upper bound for the net 
surface currents. On the other hand, the sum rule argu- 
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merit in the previous section claims that surface currents 
should be exactly canceled between opposite surfaces in 
a cylinder. 

The above observation is applicable either with or 
without the symmetry breaking fields. In the following, 
we discuss some subtleties related to the external fields 
and the thermodynamic limit, which have also been over¬ 
looked in the previous studies. Here, we focus on possi¬ 
bilities of total currents of order in d-dimensions), 

which are bulk currents but not surface currents. Firstly, 
we consider a case where H itself does not have time- 
reversal symmetry without symmetry breaking fields, al¬ 
though such a Hamiltonian would be artificial. In this 
case, we could keep L to be some large but hnite val¬ 
ues, so that the total current term becomes dominant in 
Then, we can safely apply the original Bloch- 
Bohm’s argument to conclude that total current of order 
must vanish for sufficiently large L. 

However, if we include the symmetry breaking fields 
for time reversal symmetric id, we have extra terms in 
SE^^x arising from Aiiex and the resulting 0{L) term 
in SEl x is not determined by the total current only, 
SEl,x ~ {±2Tr/L){J2[jyi + A(- • •)]) where A(- • •) repre¬ 
sents Xiijji — yj)c\/S.Cj. At this point, what is forbidden 
by the Bloch-Blohm’s argument is that the above expec¬ 
tation value is some finite value of order L/^ for large 
L. We do not know, however, which contribution in the 
bracket becomes dominant for given L, A, although the 
second term may be smaller for sufficiently small A. Fur¬ 
thermore, since we are interested in the spontaneously 
symmetry broken states, we should take the limit of van¬ 
ishing external fields A —>■ 0. When the system is defined 
on X, y = —L /2 -|- 1 , • ■ • , T /2 with an even L, we consider 


SSr] 


lim lim 

A4,0 L'I'oo 


SEl,X 

L" 


(7) 


The contribution arising from Aiiex will vanish as A —0 
due to the prefactor A in front of iiex, as long as the ex¬ 
pectation value limi-|-oo is 

not singular at A = 0. Because we have assumed that 
time reversal symmetry is broken and the corresponding 
order parameter is finite in the thermodynamic limit, this 
term should be some constant which is independent of A 
in the limit A —>■ 0, and therefore we can safely take the 
limit. It is noted that difference in the thermodynamic 
energy density 6 e 2 vanishes in this limit, while dei can be 
negative if the total surface current per volume is 0 ( 1 ) for 
finite L and is non-zero after taking the limit. However, 
we should be careful about meaning of the possible non¬ 
zero Sei- In the thermodynamic limit, quantum states 
for fixed A may be defined as 


a;o(-) = lim (Ol,a| ■ |0l,a), 

Ltoo 

(8) 

we(-) = lim {0L.x\ ■ \0 l.\) 

Ltoo 

(9) 


for local operators^. There are some subtleties in the 
thermodynamic limit, where there is no local operator 


describing the total current density for the whole system 
corresponding to fci ^ jyj)/ L'^ . In this case, the two 
thermodynamic states ojo,oJe become identical; 

uJo[A)=we[A) (10) 

for any local operator A. For example, LOg{c\iCjii) = 
AT[LL^oaE^^y^~y^\ 0 L,\\c\iCjv\ 0 L,\) = ujoichcjii) because 
lUi — Ujl < It is finite and limL-j-oo = 1. (Note 

that expectation values of U(l) breaking local operators 
vanish trivially for both states.) This is true even when 
the two states are orthogonal for finite L, and the two 
orthogonal states can converge to a single state in the 
thermodynamic limit. 

Such a behavior arises from global nature of the vari¬ 
ational state The two states |0 l_a), I^l.a), are al¬ 

most identical locally and their difference appears as a 
sum of these tiny local differences. In the present sys¬ 
tem, to obtain different states in the thermodynamic 
limit, we need to restrict the twist only for a finite sup¬ 
port such as D = {(x,y)|l < x — x,y — y < L'} where 
(x, y) is an arbitrary site. We then define a new vari¬ 
ational state \d'i^ x) — Ug\0L,\), Ug = exp[i9'yjrij] 
where O' = 2tt j L' and the summation is taken only for 
the above finite domain D^. It is noted that, in order for 
U'g to be well-defined, L' should be L/L' =integer. When 
we take the thermodynamic limit, we keep L' constant 
but increase L only. Nevertheless, we can take L' to be 
much larger than the hnite hopping range of the model 
It <C L'. Then, the energy difference for 0' = E27tIL' is 
SEl,x = ±{27r/L') E'(0L,A|jy* + A(- • • )\0l,x) + Oik/L'). 
By taking the limit, the energy difference for 0 = ±2Tr jL' 
becomes 

27r ' 

lim lim SEl^x = ±77 hm lim (0 l,a| ^ + 0{L'°), 

A4,0 L'I'OO L/' A4,0 L'I'oo * ^ 

i 

( 11 ) 

where the 0{L"^) = 0 ( 1 ) term come from the cos-term. 
Note that the hrst term does not describe total current of 
the system, but it corresponds to current running within 
the hnite region D. Here, we recall that the ground 
state of Hx in the thermodynamic limit is dehned so that 
wo(At[FlA,A]) > 0 is satished for any local operator A. 
If we take A= U'g, this means 

limu!o(Ug HxUa — Hx) = lim lim 6 E[, a > 0. ( 12 ) 

A4,0 AiO Lfoo 

By comparing Eqs. (ED and ED, we conclude that 
the current running in the domain D cannot be of order 
0{L'‘^) for sufficiently large L' ^ It- This statement is 
a bit stronger than the original Bloch-Bohm’s argument, 
since the position of D characterized by the site (x, y) can 
be arbitrary in the inhnite system. Hence, we conclude 
that there is no macroscopic flow anywhere in a thermo¬ 
dynamic system. We will see that this statement holds 
for U(l) symmetry broken systems as well in Sec. lHIl The 
only remaining possibilities would be locally circulating 
currents in an atomic scale and surface currents which 
are of smaller order in L' as discussed before. 
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C. numerical simulation 

In order to confirm the sum rule argument, we per¬ 
form numerical calculations of surface currents in a sim¬ 
ple toy model defined on a square lattice cylinder. We 
impose the open (periodic) boundary condition for the 
a;(t/)-direction. As typical examples of surface pertur¬ 
bations, we investigate effects of surface roughness and 
surface potentials. We consider the following Hamilto¬ 
nian 

Har = Hxm + HsuR, (13) 

Haro = + i42Cj2 -b + (h.c.)] 

NN 

+ '^M[n^i -ni2], 

where = ±1, and Vu is the sur¬ 

face perturbation. This model is a spinless version of 
the Bernevig-Hughes-Zhang model, and does not have 
time reversal symmetry^i^. In the present calculations, 
Vu is finite only at the left surface sites i = {x = l,y). 
We study two cases: (i) take real random values in 
[—in case of surface roughness, and (ii) I^i _2 are 
constant, {Vii,Vi 2 ) = (Vb,0), as a particular realization 
of surface potentials. Since the t'-term is simply a part 
of kinetic term which arises from spin-orbit interaction in 
the original Bernevig-Hughes-Zhang model, the current 
density operator for the /r = x, y-direction in the present 
model is given by 

jfii — — c|'_|_^j^Cii] + it[c|2Ci+/i2 ~ c|_|_^2Ci2] 

-b [cliCi+fi 2 + + (h.c.). (14) 

For simplicity, we fix (t', M) = (0.2t, —3t) as an example, 
and filling is n = ni -b n 2 = 0.7 for which the system is in 
a metallic anomalous Hall state. Temperature is fixed at 
T = 0. System size is x Ny = 80 x 20 for the case (i), 
and numerical results are checked for other systems sizes. 
For the case (ii), with use of Fourier transformation for 
the y-direction, larger system sizes are examined. 

In Fig. [21 we show the current densities jyi for the 
random potential and constant potential at a large Vq = 
2 t^t' = 0.2t together with jyi for Vb = 0 (which we de¬ 
note jyi hereafter). We see that jxi vanishes everywhere 
in the system and ly is unchanged by the left surface po¬ 
tentials. For the random potential, we take a disorder av¬ 
erage (jyi)av and then take an average of them along the 
y-direction, {1/Ny) It is seen that, although 

there are some small oscillations in the bulk region due to 
metallicity in the present model, jyi is localized around 
the surfaces. The current densities are modified from 
A around the left surface x = 0. However, we find that, 
both for the random potential and constant potential, the 
left surface current ly = (1/JVy) X^iesurf Jy* unchanged 



FIG. 2: Current density jy near the left surface for Fo = 0 
(red), random potential with Vo = 2t (green), and constant 
potential with Vo = 2t. 
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FIG. 3: Normalized histogram of ly for disorder potentials 
Fo = 0.5t (left panel) and Vo = 2t (right panel). 

from 1° = {l/Ny) E,esurf iyi. which is confirmed for sev¬ 
eral different system sizes. Conservation of the surface 
currents were also seen in the previous study for a con¬ 
stant surface potential in an insulating statei^. In the 
case of disorder potential, distribution of ly for different 
configurations of the potential is well localized around its 
mean value ly as shown in Fig.jSl It is noted that the dis¬ 
tribution of ly gets broader if we introduce inter-orbital 

random potentials X][^ii 2 cliCi 2 + (h.c.)] in addition to 
the intra-orbital potential (not shown). Although finite 
size effects become rather large in this case, the average 
surface current is still unchanged by the surface disorder. 
We have performed similar calculations for other realiza¬ 
tions of surface potentials, and confirmed that the surface 
current is unchanged by them. These numerical calcula¬ 
tions indeed support the sum rule argument in Sec. Ill Al 


D. Revisit of Orbital Magnetization Formula 

As mentioned in Sec. Ill Al the sum rule is helpful not 
only for basic understanding of surface currents but also 
for practical calculations of them. Here, based on an 
observation of the sum rule, we rederive the formula of 
orbital magnetization M at T = 02r— for a thermo¬ 
dynamically large but finite size system. We consider 
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FIG. 4: Schematic picture of the system with the confine¬ 
ment potential V which is zero in the shaded region and infi¬ 
nite in the white region. V changes in a length scale which 
is much longer than the lattice constant a. These length scales 
satisfy a <C <S i <S T. 


a non-interacting model of size L x L with the peri¬ 
odic boundary condition, whose Hamiltonian is expressed 
in a Wannier function basis as Hq = <PICiIj = 

J2ijw■ The chemical potential has been in¬ 
cluded in Hq. Then, we introduce a confinement poten¬ 
tial, Hsurf = which smoothly varies in space and 

is zero inside a region I x I while infinitely large outside. 
The length scales can be taken as 

a ^ ^ / C T (15) 


where a is the lattice constant and characterizes the 
spatial variation length of V. The system configuration is 
schematically shown in Fig. 01 Although such a confine¬ 
ment potential would not be a realistic surface potential, 
the sum rule guarantees that the surface current is equiv¬ 
alent to that with a realistic confinement potential. It is 
noted that, as long as the length scale of the confinement 
potential is much shorter than the system size ^y <g; I, 
the locally circulating current around each atom in the 
bulk is not affected by the potential. 

Generally, the paramagnetic current density in unit 
of e and orbital magnetization in a hnite system with 
boundaries are given by— , 


j{r) = 


Mz X vol = 


i / (fx[rxj(r)]^ 

^ J finite 

d^x{r-R,) xj{r)\^ 


% 


(16) 


(17) 


where Vi is a unit cell with its center position iiy The 
first term in arises from locally circulating current, 
while the second term is due to the surface current. 

We first consider the surface contribution. In the pres¬ 
ence of the smooth confinement potential, derivative ex¬ 
pansion in the Wigner representation would be legiti- 
matei^i^, where the site index i can be considered as 


a continuum variable in the length scale ^y a. The 
Green’s function is approximated up to the lowest order 
with respect to derivative of the confinement potential 
by, 


G[X,k) = GQ 


— Gq 


. 2 dXf, dk^ 

r 

. 2 dx^ dk^ 


1 dG-^ dGp - 

2 dkfj, dX^. 

2 0 dkf, dX^ 

(18) 


where Gq{X, k) is the matrix inverse of [iuj — Hp^k) — 
H(X)] = Gp^ with respect to the indices 1,1'. X = 
{xi -I- X 2)/2 is the center of mass coordinate and fc is a 
wavevector corresponding to the relative coordinate Xi — 
X 2 . In the above, we have used 0 = d{Gp^Gp)/dk^ = 
{dGp^/dkfi)Gp + Gp^{dGp/dkfi). The surface current 
along the x-direction is given by 


4 


r 

lx}. 




'"•Tic*'*" 


■ OGp^ dGp 

. dkx dky 


dGp^dGp i 
dky dkx J ’ 
(19) 


where dVjdXx = 0 near [010] surface and / dXy 
is restricted around the surface where V{Xy) = 0 
and V{Xy) = oo. tr represents summation over all 
the indices other than X. By using Gp[X,k) = 
J2n \ukn){ukn\/{iixi - Efcn - H(X)), we Can perform in¬ 
tegral over uj and simplify the expression as was done in 
Chen and Lee— to obtain 


4 


n 


z 

kn 


^ kn 

1 

/ def{e)ni^ 

kn 


1 

Z 2 


£fen<0 


z 

kn ’ 






dul^ dUkn 


1 1^2 dkfj. dky 


( 20 ) 

( 21 ) 


where / is the Fermi distribution function at T = 0. This 
is nothing but the surface current evaluated within the 
quasi-classical wavepacket theory in the previous stud¬ 
ies — The orbital magnetization arising from the 
surface current is, 

= ^ E ( 22 ) 

ekn<0 


Next, we consider the bulk contribution, = 

E,(l/2/^) /„, d^x[{ry-R,) xjl ^ W(1/2Z2) d4[r, x 

j]z where vp is the unit cell with its center Rp = 0 and 
JVi is the number of unit cells inside the confinement 
















































potential. P = vqNi is satisfied. Here, the position r 
has been replaced by = sm{qrf^)/q where q = 27r/L 
which is consistent with the periodic boundary condition. 
However, in order to calculate this contribution, we can 
safely approximate as ~ r + 0{vq /L) in the inte¬ 
grand. Besides, we can simply neglect effects of the con¬ 
finement potential and approximate the Green’s function 
as G{xi,X 2 ) — Go{xi,X 2 ) for Xi,X 2 G fo. Then, 


Therefore, we obtain 


^bulk ^ 


2 vo 


s y 




dku 


dK 


( 27 ) 


which agrees with the previous work s By collecting 
the two contributions, and we finally end 

up with the OM formula 


- — f (fx[r X {Vi -V 2 )Go{xi,X 2 )\xi=X 2 ]z 

4tovo 

(23) 


can now be directly calculated, e.g. by the 
Green’s function method in the first principles calcu- 
lationsi. Alternatively, we can also use Gq{xi,X 2 ) = 
Y.kn^kn{xi)(j)l^{x 2 )/{iuj - Ekn) for Xi,X 2 G Vq and ex¬ 
press this contribution in terms of Bloch functions by 
a formal calculation. By denoting p^, = —id^/m and 
Hk = we have 


J Vo 




in 








in 


J Vq 

/ cfx(l)l,^,p^(j)kn 
^-Skk'Snn'idk'^ f (Pxul,^,e~^’^' 

n' “'’'0 

[ Sxe-^^^idk,uln{x) 

J 7;o 

[^ ^ (j)k'n'{x)(j)l,^,{x')p^(j)kn{x') 

^ k'n^ 

/ OH 

dPxidk^ul^-^Ukn- (24) 


f j 2 9ul dHk 

IQ Ok^ Oku 

f j 2 9ul^ dSkn 

J Okfx Oku 

above expression vanishes 


(25) 

after 


j2 

d^X Ukn 

Okjj, Oku 


f f q2dul,^^ dSkn 
Jvo dk'^ Oku 


■kn 


zn jf 

IG . ■ 


M, = -k 




Ghn<0 


cPk 


Oi I * 

d\^{Hk+eku) 


dukn 


We note that, since is independent of the chemi¬ 

cal potential in insulators, we can reproduce the Streda 
formula in terms of the surface current onlyS*^, 


_ 1 

dp dp 27r^’ 



(29) 

(30) 


where p is in the gap. 

Finally, let us briefly discuss relations of the present 
results in bounded systems to the previous calculations in 
periodic systems— d^d^d^ . In the following, for simplicity, 
we consider electromagnetic coupling up to the first order 
in B = (0,0, H) at a fixed gauge, in order to discuss 
OM at zero held. For a uniform magnetic held A = 
{\/2)B X r in a bounded hnite size system, 

Hem = -[ d^xj ■ A =-(I- f (fxrxj\-B 

J finite ^ ^ J finite ' 

(31) 

holds as an operator identity, where j is the paramag¬ 
netic current Eq. (Hl-S. Expectation value of the above 
integrand is not uniform in space for both expressions, 
and the surface current contribution is localized around 
the surface. Nevertheless, the integrated energy (Hem) is 
independent of surface conditions, since the OM is a bulk 
quantity as implied by the sum rule argument. Then, it 
is quite natural that the total energy including (Hem) 
of the bounded system is equivalent to that in a peri¬ 
odic system of the same volume with the uniform mag¬ 
netic held B. This should be true from a macroscopic 
point of view that total energy of a system which is an 
extensive quantity does not depend on boundary condi¬ 
tions in the leading order of the system size. Therefore, 
OM calculated by derivative of the free energy with re¬ 
spect to B in the periodic system is equivalent to that 
in the bounded system computed either from Eq. (II3 or 
from derivative of Eq. dSB). However, from a microscopic 
point of view, the coincidence of these two quantities is 
not a priori guaranteed, because OM in the periodic sys¬ 
tem is dehned only by the derivative of the free energy 
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but is not given by an expectation value of an OM op¬ 
erator such as Eq. dlB since it is ill-defined under peri¬ 
odic boundary conditions. In this sense, the sum rule of 
surface current density gives a microscopic basis for the 
macroscopic equivalence of total energies under different 
boundary conditions. 

III. SYSTEM WITHOUT U(l) SYMMETRY 
A. general argument 

In this section, we consider superfluids without time 
reversal symmetry. As noted briefly in the previous sec¬ 
tion, in case of systems without U(l) symmetry, we can¬ 
not apply sum rule arguments on robustness of surface 
currents against surface perturbations. The main diffi¬ 
culty arises from the modification of the velocity matrix 
in the Green’s function formalism. 


matrix such as spin-orbit coupling and inter-orbital hy¬ 
bridization. Similarly to the U(l) broken case, we will 
have the same problem and cannot apply sum rule argu¬ 
ments to those cases. 

Similar difficulty arises in the Bloch-Bohm’s argument 
for superfluidity. Although it is not helpful for sur¬ 
face currents as was discussed in Sec. iim we briefly 
discuss it in superfluids for a comparison with U(l) 
symmetric systems, focusing on macroscopic bulk cur¬ 
rents. The largest difference between U(l) symmetric 
systems and U(l) broken systems comes from the exter¬ 
nal field Ai/ex under twist by Ue = exp[i0 ^ UjiT-j]. When 
XHex = + (h-C-)j if is transformed as 

ulXH^^Ue = + (h.c.). (34) 

If we evaluate energy difference between the ground state 
| 0 l,a) and a variational state | 0 l,a) = Ug\0L,\), we ob¬ 
tain, in the leading order oi 0 = 27rn/T, 


Vy{ky) Vy{ky) = 


1 


Vv{ky) 

0 


y 2 ^ 

= -Q^ 
2^ dk,. 


-Vyi-ky) 


(32) 


where Q = diag(l,—1) and K{ky) = 
d:mg{k{ky),-k^{-ky)) in the Nambu space. The 
charge matrix Q is not a unit matrix, since the charge 
carried by the particles is assigned to be - 1-1 while it is 
— 1 for the holes. Correspondingly, is not simply 
given by a simple form as in Eq. Q. For this case, we 
cannot simply perform the summation over ky in Iy°^, 
and cannot obtain a sum rule for the surface current. 
Instead, we can consider a related quantity 1*°* which is 
given by 


^y°^ = [*5 ^v{ky)G{k) 


N, 




dG 


-1 


dk„ 


-G 


1 V- 1~ 
^ ^dky . 


• (33) 


Here, G and E are respectively the Green’s function and 
the selfenergy in the Nambu representation, and trace in¬ 
cludes summation over the Nambu space. While de¬ 
scribes difference between particle-like contribution and 
hole-like contribution, is related to an equal-weighted 
sum of these contributions and is conserved in the pres¬ 
ence of surface perturbations due to the sum rule. This is 
analogous to the problem of Fermi surface volumes where 
each of them and hence difference among them are not 
conserved in general, while their total sum is unchanged 
by interactions as stated in Luttinger’s theorem. We note 
that, for general S\J{N) currents such as spin currents 
and orbital currents, corresponding charge matrices are 
not the unit matrix in the spin/orbital space and there 
would be off diagonal matrix elements in the velocity 


SEl,x = {6l,\\Hx\6l,x) — ( 0 l,a|Ha| 0 l,a) 

= A^(cOs[6l(?/i -bj/j)] - l)(0L,A|cJ;Ai;j7'c];,|0L,A) 

-I- A ^ f sm[ 6 »(y, -b yj)] {0L,\\4Aaji’cli> |0l,a) 

-b (h.c.). (35) 

Since yi + yj can be of order L, we cannot Tayler ex¬ 
pand cos[' • ■ ]/ sin[' • ■ ] and neglect higher order terms in 
(yi + yj)/L. On the contrary, above two terms will be of 
order and because (cos[ 0 (?/i-b?/j)] — 1 ) = —2siT?[9{yi + 
yj)/2] < 0 while sin[ 0 (yi -b yj)] is oscillating in sign, the 
first term would become dominant. Indeed, the latter 
term will vanish if |0l,a) is translationlly invariant in a 
long distance scale ^1/0. Sign of the first term can be 
evaluated, once we simply assume ( 0 L,A|Ai?ex| 0 L,A) < 0 , 
which is reasonable for spontaneous symmetry breaking. 
This assumption is a variant of the statement that ex¬ 
ternal magnetic fields parallel to the magnetic moment 
lower the total energy in conventional ferromagnets. In¬ 
deed, ( 0 L,A|AI?ex| 0 L,A) is & part of condensation energy 
of superfluidity, and therefore should be negative when 
the finite size system has instability towards the corre¬ 
sponding superfluidity. If the above assumption really 
holds, we see that SE^^x ^ X x o{k) > 0 in the leading 
order of L by noting that —2sik[9{yi + yj)] < 0 can be 
approximated by a negative constant of order unity. 

However, similarly to U(l) symmetric systems, the 
variational state ujg becomes identical to ujq in the 
limit A ^ 0, and Ymix\a^^^L^oo ^El^x /vanishes be¬ 
cause of the prefactor A in front of i?ex- Therefore, 
we introduce the other variational state ]0']^ x) which is 
twisted only in a finite domain D. By repeating the 
same calculation, we find that the leading contribution 
in SEl^x comes from AHex which is 0{L'^), and next 
leading 0 (L')-contribution is the total current term if 
^ 0{L'‘^). For the twist only in D, how¬ 
ever, the former 0{L'^) term will vanish in the limit 
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liniA 4 .o liHiL-t-oo and we obtain Eq. (ITT|) as in U(l) 

symmetric systems. (As mentioned for U(l) symmetric 
systems, limi.|.oo would not be sin¬ 

gular at A = 0 and will vanish for A —?► 0.)Therefore, we 
arrive at the same statement as for U(l) symmetric sys¬ 
tems that macroscopic currents are not allowed anywhere 
in the ground states of superfluids iir the thermodynamic 
limit. This statement also holds in the presence of static 
magnetic fields in superconductors and macroscopic su¬ 
percurrents flowing in the bulk are not allowed at equilib¬ 
rium. Therefore, if a Fulde-Ferrell state or a helical state 
with non-zero center of mass momenta of Cooper pairs is 
realized, there should be some counter-propagating cur¬ 
rents which compensate the macroscopic supercurrents. 
For example, it was pointed out that in the helical states 
in noncentrosymmetric superconductors under magnetic 
fields, supercurrents are canceled by magnetization cur¬ 
rents and there are no currents in the thermodynamic 
limit^SilL. 

The Bloch-Bohm’s argument can predict vanishing 
macroscopic currents which is proportional to domain 
volumes, while the Green’s function approach is less help¬ 
ful for U(l) broken systems. However, neither of them 
can exclude possible net currents due to incomplete can¬ 
cellations of surface currents. 


B. numerical simulation 

Although the sum rule discussions cannot be applied 
to superfluids, it is still possible that the surface mass 
current is robust against surface perturbations by some 
other reasons. In order to investigate this, we examine 
numerically surface currents in two simple models, a non- 
chiral p-wave superfluid based on the model (1131) and a 
chiral p-wave superfluid. We focus on neutral fermions 
and do not consider Meissner effects in the present study. 
Temperature is fixed at T = 0. 

1 . non-chiral p-wave superfluid 

Firstly, we consider a non-chiral p-wave superfluid 
based on the model (USD in order to examine how U(l) 
symmetry breaking modifies the previous results in Sec. 
Ill Cl The Hamiltonian is 

I^AHSF = -SIaho + Hsnif — g’^^nurii+yu (36) 

.Hsurf — ^ ^ ^iW ^ 

where g is an attractive interaction for py-wave superflu¬ 
idity. As in Sec. Ill Cl we again consider two particular 
examples of surface perturbations, a random potential 
and a constant potential. In case of disorder potential, 
inter-orbital surface potentials Vii 2 = are introduced 
in addition to the intra-orbital potentials 1 ^ 11,22 = 

This Hamiltonian is one of the simplest models to discuss 



FIG. 5: The surface current ly normalized by Iy{g — 0) in 
the absence of surface disorder, Vb = 0. 


intra-orbital superfluidity in the model (1131) . and this su¬ 
perfluidity itself does not break time-reversal symmetry. 
We perform mean field calculations of the superfluidity. 
The mean field calculations are performed for each disor¬ 
der configuration in the case of disorder potential, which 
is repeated until averaged physical values become con¬ 
verged. It is noted that, in the cylinder geometry where 
the periodic boundary condition is imposed for the y- 
direction, there is no zero-energy Andreev bound state at 
the surfaces, which makes numerical calculations rather 
stable. The system size mostly used for the random po¬ 
tential is Nx X Ny = 40 X 20 and results are qualitatively 
unchanged for other sizes up to x Ny = 60 x 20. For 
the constant potential, similarly to the previous section, 
we can perform Fourier transformation for the y-direction 
and study larger sizes. 

We show the surface current as a function of 5 at Vb = 0 
in Fig. [5] The surface current is suppressed by the Py- 
wave superfluidity. It is noted that similar behaviors are 
also seen for spatially uniform gap functions whose ampli¬ 
tudes are chosen to be consistent with the self consistent 
calculations. For non-self-consistent gap functions, we 
can tune the gap amplitudes A 1.2 for each orbital inde¬ 
pendently in order to investigate A 1 2 -dependence of the 
surface current. By calculating the surface currents for 
such Ai _2 we see that is determined by detailed bal¬ 
ance between the gap amplitudes for the two orbitals (not 
shown). In the self consistent calculations, ratio between 
Ai 2 is determined by the gap equation, which then leads 
to the non-monotonic behavior of Iy{g) (Fig. [5]). 

In the presence of the surface perturbation potential 
Viip, the current density is modified as in the previous 
section. We find that the resulting left surface current 
ly can also be changed from 7° in contrast to U(l) sym¬ 
metric systems, while is unchanged. For the random 
potential, the surface current ly is suppressed as shown 
in Fig. [6l ly is almost unchanged up tp Vb — t, and it 
decreases by further increasing Vb. The reduction of ly 
by the surface roughness well agrees with our naive ex¬ 
pectation that disorder would generally suppress surface 
currents. On the other hand, for the constant potential 
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FIG. 6: The surface current ly normalized by Iy{Vo — 0) 
when g = 5t. 



FIG. 7: The current density jy near the left surface for the 
constant surface potential {Vii,Vi 2 ) = (Vb,0). 


along the left surface sites, 1 ^ 2 , 14 i 2 ) = (Vb, 0, 0), the 
surface current ly shows non-monotonic Vo-dependence. 
ly is quickly suppressed as Vq is introduced, and then 
it turns to increase exceeding when Vo is sufficiently 
large. In order to understand this behavior, we show the 
current density jyi in Fig. [T] As Vq is increased, jyi gets 
suppressed at the left surface sites i = {x = 0,y), but 
at the same time, it is increased at the next surface sites 
i = (x = \,y). The reduction at x = 0 sites determines 
ly for small Vq <^t, while ly is dominated by the contri¬ 
bution from X = 1 sites for large Vq. Since we now do not 
have U(l) symmetry and an associated sum rule, these 
changes in jy do not necessarily cancel out and indeed 
they add up to give non-trivial finite values. Therefore, 
ly deviates from ly and shows the non-monotonic be¬ 
havior in the present system. 


S. chiral p-wave superfluid 

As a second example of superfluids without time rever¬ 
sal symmetry, we study a chiral p-wave superfluid on a 
square lattice in which surface current is generated by the 
superfluidity itself. The problem of spontaneous surface 
current and corresponding 0AM, often referred to as “in¬ 
trinsic angular momentum paradox”, has been discussed 



FIG. 8: The chemical potential p, amplitude of the gap 
functions in the bulk A, and spectrum gap 5E for filling n = 
0.2. The purple dotted line represent pc = —4t. 


for more than 40 years— Although most of the previ¬ 
ous studies focus only on the weak coupling BCS region, 
here we discuss both the BCS region and the BEC re¬ 
gion on an equal footing. In the present study, similarly 
to the previous models, surface roughness is introduced 
as a typical example of surface perturbations. We con¬ 
sider the following Hamiltonian 

HpSF — ^ ^ 9 ^ ^ -f Hgurf, 

NN NN 

-^surf ^ ^ 

where Vi is finite only at the left surface, i = (x = l,y). 
The hopping and interaction are allowed only for the 
nearest neighbor sites, and the same cylinder geometry 
as in the previous sections is used. The system size is 
Nx X Ny = 60 X 20, and we have confirmed finite size 
effects are negligibly small for this size by performing 
similar calculations for other sizes. We perform mean 
field calculations of the superfluidity, which is a good ap¬ 
proximation even for a large g at zero temperature since 
there are no thermal fluctuations in the ground states 
The current density is simply 

dui ~ ~ ^l+p(T^i<r]- (38) 

(7 

Before discussing surface disorder effects, we first ex¬ 
amine basic properties of the model in the absence of 
surface disorder. The present model exhibits a quantum 
phase transition— when p = pc = —4t. For |^| < \pc\, 
the system is in the BCS region where there are gapless 
chiral edge modes at the surfaces, while for |/r| > \pc\, the 
system is in the BEC region where there is a spectrum 
gap. In Fig. [51 we show the chemical potential p, ampli¬ 
tudes of the gap functions at a center site of the system 
A, and spectrum gap SE for fixed filling n = 0.2. The 
quantum phase transition takes place around 5 ~ 8 t — 9t 
for this filling where p crosses pc, and we have confirmed 
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FIG. 9: The current density jy near the left surface for 
Vb = 0 at n = 0.2 with g = 5t (BCS state, green) and g = 15t 
(BEG region, blue). 


FIG. 11: The current density jy near the left surface for 
Vo = 8t at n = 0.2 with g = bt (BGS state, green) and 
g = 15t (BEG region, blue). 
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FIG. 10: The surface current for different filling n at g = 5t 
and g = 15t in unit of the effective mass mefr = l/(2ta^). 


that similar behaviors are seen for other low filling. It is 
noted that, when filling is high, n ~ 1, the ground state 
stays in the BCS region even for large g and it is hard to 
realize the BCS-BEC phase transition. 

We show the current density jy at u = 0.2 for g = 5t 
(BCS region) and g = 15t (BEC region) in Fig. |9] as an 
example, jx vanishes everywhere in the system. The 
current density in the BEC region is more strongly lo¬ 
calized near the surface than that iir the BCS regioir, 
and it oscillates in sign dependiirg on the distance from 
the surface. As a result, the surface current in the BEC 
region is smaller than that in the BCS region in the 
present lattice model, as shown in Fig. [Tni We see 
that the overall behavior of in the BCS region as a 
function of filling n is consistent with the recent worb2£. 
The surface current for low filling and weak coupling 
limit approaches I = n/(4meff) where rueff is an effective 
mass TTieS = l/{2ta^) with the lattice constant a. Un¬ 
der an assumption that the surface current is constant 
along a boundary of a finite system, this gives the 0AM 
= rUeS §[rxl]zdl = N/2 where N is the total number 
of fermions in agreement with the previous calculations 
for continuum systems without lattice potentials— 

In the present lattice model, in contrast to the contin¬ 
uum systems where = N/2 holds both in the weak 


coupling region and strong coupling region, ly and the 
corresponding 0AM is decreased when the coupling 
constant g is increased. This is because, even for low 
filling n <g; 1 where lattice effects is expected to be less 
important, the smallest size of a bosonic molecule of two 
fermions is bounded by the lattice constant for a non-s- 
wave superfluid on a lattice, and therefore, presence of a 
lattice is significant especially for the BEC region rather 
than the BCS region. Because of this lattice effect, the 
surface current per fermion for strong g = I5t is almost 
independent of filling as seen in Fig. [Tni although the 
system stays in the BCS region for high filling n ~ 1. 

Now, we discuss effects of the surface disorder. The 
current density jy averaged over disorder configurations 
for n = 0.2 at a large Vb = is shown in Fig. [TTJ For 
the BCS region, jy is suppressed especially at the surface 
sites 1 = 0 and it is enhanced at inner sites a; = 3,4 
to partly compensate the reduction, while jy in the BEC 
region becomes strongly oscillating and effects of U prop¬ 
agate into further inner sites a; = 4, 5. As in the previous 
model (1571) . change of jy does not need to obey a sum 
rule and the surface current can be modified from ly. 
Indeed, as shown in Fig. [TTl the surface current is sup¬ 
pressed by the surface disorder. The decrease of l/j in 
the weak coupling BCS limit is consistent with the pre¬ 
vious studie o^^'^^ 1 ^^ . Interestingly, l/j decreases not only 
in the BCS region but also in the BEC region. This 
would be because the smallest bosonic molecule size is 
bounded by the lattice constant and the disorder length 
scale is of the same order in the present model. We have 
confirmed similar behaviors of ly for different parameter 
sets {g,n). For a comparison, we also calculate surface 
currents for non-self-consistent gap functions which are 
constant in space and whose amplitudes are chosen to be 
consistent with the self-consistent calculations. The sur¬ 
face currents are suppressed in a similar way as in the self 
consistent calculations, which means that change of the 
gap functions around the surface by Hsmf is not impor¬ 
tant for the reduction of ly. We note that, although the 
reduction of ly by the surface disorder is moderate, it 
was pointed out that, for domain boundaries with oppo- 
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FIG. 12: The left surface current ly for different Vo at g = 
5t (green square) and g = 15f (blue circle) when filling is 
n = 0.2. Filled symbols correspond to non-self-consistent gap 
functions whose amplitudes are chosen to be consistent with 
self consistent calculations at g = 5f and g — 15t, respectively. 


site chiralities in the BCS region, the boundary current 
strongly relies on boundary conditions and it can change 
even its directional^. 


C. discussion 

We have discussed suppression of surface currents in 
cylinder systems in the previous sections. For a realis¬ 
tic finite size system with open boundary conditions for 
all directions, surface current flows along the surface and 
generates global rotation. In the absence of surface per¬ 
turbations, its magnitude is obviously uniform along the 
surface. If some parts of the surface are perturbed, the 
surface current would be changed not only at the per¬ 
turbed surface but also at the whole surface, because of 
the continuity of the current density. Therefore, associ¬ 
ated 0AM would also be modified. 

Although we have examined particular realizations of 
surface potentials among possible surface perturbations, 
existence of surface perturbations which change ly con¬ 
ceptually distinguishes a system without U(l) symmetry 
from a system with U(l) symmetry. The surface current 
is not uniquely determined as a bulk property and there 
may exist surface perturbations which drastically change 
it in the former, while the surface current is an intrin¬ 
sic quantity in the latter. The absence of a sum rule 
for the surface current density and the numerical results 
suggest that there is no bulk-surface correspondence for 
surface currents and corresponding 0AM in superfluids, 
and surface conditions should be fully taken into account 
in order to calculate these quantities in contrast to some 
of the previous studies for chiral superfluids From 

these discussions, it is considered that surface currents 
and 0AM in superfluids with broken time reversal sym¬ 
metry would be subtle quantities, and experimentally, 
one needs to control surface conditions carefully in order 
to measure these quantities. 

Let us briefly discuss the change of surface currents 


and 0AM in superfluids in view of thermodynamics. In 
the present study, we have implicitly assumed that the 
lattices (or containers for non-lattice systems) are at rest 
in the laboratory frame by some reasons. If the lattice 
or container is fixed spatially to a much larger environ¬ 
ment, the system composed only of the superfluid and 
lattice does not conserve the 0AM and it is an open 
system with respect to angular momentum. Although 
0AM of the superfluid alone can be changed by surface 
perturbations, total 0AM of the whole system includ¬ 
ing the large environment should be a conserved ther¬ 
modynamic quantity and is independent of surface de¬ 
tails of the lattice/container. If the system composed 
of a superfliud and a lattice/container is suspended in 
the midair and set to be at rest, as temperature is de¬ 
creased down to the superfluid transition temperature, 
the lattice/container should start to rotate in an oppo¬ 
site direction to the superfluid rotation in order to keep 
the total 0AM of the whole system zero due to the an¬ 
gular momentum conservation. This is analogous to the 
Einstein-de Haas effect and each of the 0AM for the su¬ 
perfluid and lattice/container would depend on surface 
conditions in the present system. 

In the present study, we have not taken into account 
the electromagnetic field which couples to charged par¬ 
ticles. Indeed, it is especially important in superconduc¬ 
tors which exhibit Meissner effect. If Meissner effect is 
included, current density distributions are modified and 
net surface currents would vanish for uniform supercon¬ 
ducting state o^^i^^ . The problem of Meissner effect would 
be more complicated in a system where there exist circu¬ 
lating currents even in non-superconducting states, such 
as the model dSZl) and ferromagnetic superconductors. 
This issue is left for a future work. 

When we were finalizing the present paper, we became 
aware of a relevant article by Kusama and Ohashi— 
which claims that, when surface currents are not canceled 
between left and right surfaces in a cylinder, supercur- 
rent will compensate this and total current will vanish 
in superfluids. Although this is an important possibil¬ 
ity for a vanishing total current, this issue has not been 
well understood. For example, in Kusama and Ohashi—, 
the current density even near the non-perturbed surface 
is strongly changed from the original configuration if su¬ 
percurrent is included. This seems unphysical, since ef¬ 
fects of local perturbations only on a surface should not 
propagate to the opposite surface. Secondly, although 
they compared the free energy density of different sys¬ 
tem sizes for a technical reason, this cannot be justified 
in general. Their discussion is motivated by the original 
Bloch-Bohm’s theorem which is not helpful for surface 
currents as shown in Sec. Ill B[ and further investigations 
would be required to understand this issue. Here, in or¬ 
der to have an insight, let us breify consider a possible 
supercurrent in a cylinder L x L where open (periodic) 
boundary condition is imposed for a;(?/)-direction, and 
surface perturbations are introduced only for one surface. 
Supercurrent density is roughly proportional to V^(a;, y) 
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for a gap function A = e®‘^A where A is a gap func¬ 
tion without a modulation. In the cylinder considered, 
the supercurrent density should be translationally 

symmetric along the y-direction. Therefore, (p must be 
(f> = QxX + qyU with = 27rn/L, [n = 0, ±1, • • •), and 
qx should be zero for a vanishing total supercurrent in 
the s-direction. Besides, qy should be of order 2tt/L so 
that the current density near the non-perturbed surface 
remains unchanged. In such a case, qy ~ 27:jL intro¬ 
duces supercurrent density ^ 1/L at every site, resulting 
in a total supercurrent ^ L in d-dimensions) in the 

whole system, which is the same order as the surface cur¬ 
rent. We note that, however, two thermodynamic states 
constructed from wavefunctions with A and A respec¬ 
tively could not be distinguished by local operators and 
therefore they converge to a single state as in Sec. IIIBI 
Besides, if we consider a semi-infinite system and discuss 
it within weak coupling approximations as in the previ¬ 
ous studies^i^, we could not impose a boundary condi¬ 
tion at infinite with an infinitesimal supercurrent density. 
Furthermore, if we consider a realistic finite size sample 
with boundaries and introduce surface roughness, it is 
impossible to realize a uniform supercurrent density and 
possible current density configuration especially near the 
surface would be quite complicated. Therefore, possible 
compensation by supercurrent is a subtle issue. In order 
to discuss such a subtle issue, we would also have to be 
careful about validity of the mean field approximations. 
Further investigations of gap functions and supercurrent 
may be required for clarifying a role of supercurrents. 


IV. SUMMARY 

In summary, we have investigated equilibrium surface 
currents in systems with or without U(I) particle num¬ 
ber conservation. For the systems with U(I) symmetry, 
we showed that the surface currents are independent of 
surface perturbations based on the sum rule for current 
densities, which was confirmed by numerical calculations 
for a concrete model with surface perturbations. There¬ 
fore, the surface currents and corresponding orbital mag¬ 
netization are bulk quantities which are robust against 
surface conditions. The sum rule argument is also ap¬ 
plicable to the Landau diamagnetism, and it would give 
a new understanding on the known equivalence between 
the bulk approaches and the surface approaches. On the 
other hand, in superfluids which do not have U(I) sym¬ 
metry, the surface currents are changed by surface per¬ 
turbations. Especially, in a chiral superfluid on a lattice, 
the surface current is suppressed by surface disorder not 
only in the weak coupling BCS region but also in the 
strong coupling BEG region. These results imply that 
surface mass currents and orbital angular momentum in 
superfluids with broken time reversal symmetry would be 
subtle quantities and depend on surface details. Experi¬ 
mentally, one needs to control surface conditions carefully 
in order to measure these quantities. 
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Appendix A. Sum Rule for Continuum Model 

Our sum rule argument holds also for continuum mod¬ 
els with lattice potentials. We consider a general Hamil¬ 
tonian defined on a cylinder L x L with the open (peri¬ 
odic) boundary condition for a;(y)-direction, 

H = [ (fxtp^lC’tlj + Hint + HsuTf, (Al) 
Jl^ 

where K, is the single-particle Hamiltonian including 
the lattice potential and tp{r) = (^|(r),'0|(r)) is the 
fermionic field operator. Then we expand the field oper¬ 
ator in terms of non-interacting single-particle wavefunc¬ 
tions. 


V>kynir)ck^n, (A2) 

kyU 

— ^kyn^kyUt (-^^) 

where ky is Bloch wavenumber along the j/-direction and 
n represents other indices including a quantum number 
corresponding to position x. The Matsubara Green’s 
function G{r,r') = is also expanded as 


Gir,r',iuj)= ipkynir)gnn'iky,iuj)<pY,{r'), (A4) 

kyun' 

ff(fcyGw,) = - S]"\ (A5) 


where (go)„rl' = “ efc„n]^nn' is the non-interacting 

Green’s function and E is the selfenergy in the (pk n-basis. 
Similarly to the lattice models, spontaneous symmetry 
breaking order parameters are easily incorporated into 
g. In the cylinder, the surface current averaged over the 
g-direction is simply given by 


rL(R) 

y 


I 

L 

I 

L 



Y [ ^^2;gy(r), 


(A6) 


where jy is the current density. As discussed in the 
main text and Ref. 1^ . contributions to the surface cur¬ 
rent come only from ^^(A/j) and not from bulk regions. 
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Therefore, when U(l) charge symmetry is present, the 
total surface current = {ly + ly) is written as 


Ttot _ _ 

y ~ L 


'L 2 


d^x{jy{r)) 


L Jj ^2 2m y 


[dy> - dy]G{r,r',T = 0 )|r.=r 




9go \k) 

dky 


9{k) 


(A7) 


Here, trace describes summation over n and w. We can 
now follow the same argument as in the main text, and 
show = 0 even in the presence of left surface per¬ 
turbations. This means that the surface current is un¬ 
changed by surface perturbations. 

It is noted that the surface currents in lattice models 
are obtained by the following replacement in Eq. (IA6I) . 


where wpi is a Wannier function and (• • • = f (Px. 

This replacement is verified when the chosen Wannier 
function is well localized in a length scale which is much 
smaller than the system size. By similar replacements, 
other local site quantities such as particle density at 
site i become equivalent to the usual Wannier basis de¬ 
scriptions, 

= T,iCR,iCRii- Even 
when the replacement of the surface current operator is 
legitimate, however, the resulting surface current might 
depend on Wannier functions or gauge of Bloch func- 
tions^ii^, although the original definition (IA6I) is indepen¬ 
dent of them. 


f (fxij^-idy)ij 

zGsurf 

= Y 'll {wRi\- idy\wRn')y^cl^iCRn' 

zGsurf RR' .11' 

RR'Gsurf W 
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